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We describe a novel simulation method that eliminates the slowing-down problem in the Monte 
Carlo simulations of imaginary-time path integrals near the continuum limit. This method combines 
a stochastic blocking procedure with the multigrid method to rapidly accelerate the sampling of 
paths in a quantum Monte Carlo simulation, making its dynamics more ergodic. The effectiveness 
and efficiency of this method are demonstrated for several one-dimensional quantum systems and 
compared to other standard and accelerated methods. 
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I. INTRODUCTION 

The path integral formalism Q, Q can be used to de- 
scribe the statistical mechanics of a quantum system. In 
the path integral representation, every quantum particle 
maps onto a cyclic gaussian string. In most path integral 
simulations, the quantum strings are first discretized and 
then the sum over all paths is carried out by Monte Carlo 
sampling. In the discretized form, the path of a quantum 
particle is isomorphic to a classical gaussian ring polymer 
yf , and the statistical weight associated with each parti- 
cle (in the canonical ensemble) takes the following form: 



W 



exp 



2A 



(1) 



where Xj is the position of the j-th bead on the ring, 
A = etf/m, e = f3/P = (PkgT)- 1 , m is the mass, T 
is the temperature and the bead indices j are cyclic (i.e. 
P + l = l). 

To obtain the correct quantum limit, the continuum 
limit P — > oo must be taken. In this limit, however, the 
harmonic bonds between successive beads on the ring 
polymer become very stiff. This causes the simulation to 
slow down dramatically on approach to the continuum 
limit, in a way that is very similar to a system undergo- 
ing a second-order phase transition. As a result, a Monte 
Carlo algorithm employing only local updates will equi- 
librate extremely slowly. 

This slowing-down problem in path-integral simula- 
tions can be remedied in two ways. First, one can use 
a more accurate approximation for the short-time prop- 
agator in the path integral instead of the "primitive" ap- 
proximation in Eqn.Q, with the hope that not too many 
beads will suffice to accurately approximate the path in- 
tegral Q . The second way is to devise Monte Carlo meth- 
ods which employs nonlocal updates to hopefully remove 
or at least ameliorate the slowing-down problem. 

In this paper, we are concerned with the second ap- 
proach. Several alternatives to the canonical single- 
particle Metropolis Monte Carlo method p| have been 
proposed for this purpose. First, there is the so-called 



"staging" method ||J [D, lij li| • It attempts to reduce the 
correlation among the beads on the polymer by trans- 
forming to a new set of coordinates which diagonalizcs 
the kinetic part of the action. However, in the presence 
of a nonzero potential V, the transformed coordinates be- 
come correlated to each other through V again. To trun- 
cate these correlations, the staging method performs this 
transformation for one short segment of the ring at a time 
(hence the name "staging"). Updates in the transformed 
coordinates are done by direct sampling from indepen- 
dent distributions, but the new coordinates are accepted 
or rejected together based on a Metropolis criterion for 
the potential part of the action. A second method, con- 
ceptually similar to the first one, is known as the Fourier 
path integral method [13, [H 111 El Here the Monte 
Carlo moves are performed in the Fourier modes of the 
path, which also diagonalizcs the kinetic part of the ac- 
tion. But the Fourier modes are also correlated through 
V. A third method, which is also based on similar ideas 
as the first two, is the bisection methodf!, Il4|. Instead 
of generating the Fourier modes of the path, the bisec- 
tion method samples the midpoint of a large segment 
of a free-particle path and accept or reject it using an 
approximation for the long-time action that includes ef- 
fects of the potential. If the midpoint is accepted, the 
midpoints of the two shorter subscgmcnts on each side 
of the midpoint are then generated, and this process is 
iterated until every bead on the entire segment is gen- 
erated. Finally, there is a multigrid-based method, first 
applied to path integrals by Janke and Sauer 
They propose moving whole blocks of neighboring beads 
on multiple length scales using a Metropolis algorithm, 
and they cycle through the different length scales in a 
systematic way. 

In this paper, we describe another Monte Carlo 
method. This method is actually related to the four 
methods described above, but as we will show, its formu- 
lation is more general, its applications are more powerful 
and its efficiency is much higher than the previous meth- 
ods. Our method combines a stochastic blocking proce- 
dure, often referred to as the Swendson-Wang method 
, with multigrid ideas 0, in an attempt to for- 
mulate a set of equilibrium stochastic dynamics that is 
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highly ergodic for path integral simulations. The idea 
for this type of multigrid method was first proposed by 
Kandel et al. for an Ising model at criticality 20] . 

The method will be described in Sections ITT1 IIIII Sec- 
tion [H] provides the general concept of the method, and 
Section ITTT1 applies the concept to the path integral prob- 
lem. Section IIVI will compare the method against oth- 
ers for several examples of 1-dimensional systems with 
single- and double- well potentials. 

II. GENERAL FORMULATION 

The general idea of the multigrid Monte Carlo method 
has been described in detail by Kandel et al. [2£j for the 
Ising model. We will not attempt to reproduce all the 
details here. Instead, we will summarize the essentials 
in this section, using a language which is closer to path 
integrals. The specific application of these ideas to path 
integral simulations will be described in detail in Sect. IIIII 

Consider a system with partition function 




where S = ^2 a u a , and each u a is a real- valued inter- 
action term involving any number of the N particles in 
the system. Of course, all classical systems, as well as 
quantum systems that can be mapped onto isomorphic 
classical polymeric systems through the path integral for- 
malism, have partition functions of this form. The cor- 
relations among the particles arise from the interactions 
u a . 

The multigrid Monte Carlo method is based on a com- 
bination of the stochastic blocking and multigrid ideas. 
We will first describe the stochastic blocking procedure, 
often referred to as the "unigrid" method. To acceler- 
ate the dynamics of the system, the unigrid method pro- 
ceeds in two stages. First, with the current configuration 
X = {x±, • ■ ■ xn}, we attempt to remove some of the 
correlations from among the particles by "killing" the in- 
teraction terms u a one by one: For each u a , we consider 
either "deleting" it entirely from the action S with prob- 
ability pd — c a exp(u a ) or "freezing" it with probability 
Pf = 1 — pd- If an interaction is deleted, the ensuring 
simulation can update X — > X' without any regard to 
u a . If on the other hand u a is frozen, the ensuring simu- 
lation must not change the value of u a during any update 
X — > X'. To ensure that pd and pf G [0, 1], the coeffi- 
cient c a must be chosen to be smaller than exp(— it*), 
where u* a is the largest possible value for u a . After all 
the interactions have been killed (deleted or frozen), the 
particles can be divided into separate clusters - particles 
in the same cluster are connected by frozen bonds, while 
particles in different clusters are no longer correlated with 
each other. 

In the second stage of the simulation, we can update 
each cluster separately with a Monte Carlo move that 
preserves the frozen bonds inside that cluster. After all 



the clusters have been updated, we can restore the inter- 
action terms and repeat the procedure starting from the 
first stage again, or we can use a few local Metropolis 
moves to update the system before starting the stochas- 
tic blocking procedure again. It can easily be shown that 
this two-stage procedure satisfies detailed balance and 
therefore produces the correct statistical sampling^]- 

Under the unigrid method, interactions that are strong 
will more likely be frozen and those that are weak will 
more likely be deleted. This operation aims to remove 
some of the interactions from the system and this can 
potentially make the subsequent updates more ergodic. 
Whether this is actually the case will depend on two fac- 
tors: (1) whether the length scale of the unconnected 
clusters are actually small enough, and (2) whether there 
exists an efficient way to update each cluster without 
disturbing the frozen interactions. In reality, the length 
scale of unconnected clusters resulting from the stochas- 
tic blocking procedure can still be quite large, and hence 
a lot of the correlations remain in the system. In ad- 
dition, for systems with continuous coordinates, finding 
an efficient way to update all the particles in an uncon- 
nected cluster while preserving the frozen interactions 
is not always trivial. Therefore, the stochastic blocking 
procedure ends up not being as useful as it may appear. 

To repair this and to completely remove the residual 
correlations, we incorporate multigrid ideas and try to 
force the clusters to break up into smaller pieces of vary- 
ing length scales. To achieve this, the multigrid method 
first divides the particles into sets, each having a dif- 
ferent length scale. In this context, the definition of 
the "length scale" should be based on an intuitive un- 
derstanding of the physical origin of the correlations in 
the system. (For example, in path integrals, the dom- 
inant correlations among the beads on a ring originate 
from the harmonic bonds, so these correlations can be 
decomposed into a hierarchy of gaussian fluctuations on 
different length scales.) After a definition of these sets 
is made, we proceed as before but with the stochastic 
blocking procedure applied to only particles belonging to 
a single length scale. As such, we kill all the interactions 
between particles of that length scale, while the other in- 
teractions among particles on all other length scales are 
kept alive. The result of the blocking procedure produces 
clusters that are unconnected by interactions on the cur- 
rent length scale, but these "unconnected" clusters are 
not totally independent because they are still correlated 
with each other through the interactions that are kept 
alive. We update the unconnected clusters as before, but 
to maintain detailed balance, each update will also have 
to be accepted or rejected using a Metropolis criterion 
based on all the live interactions that are linked to that 
cluster. After updating all the "unconnected" clusters on 
one length scale, we can proceed to another one. In this 
manner, the multigrid method systematically breaks up 
all the remaining correlations on every length scale and 
the slowing-down problem can be completely eliminated. 
It can also be shown that this multigrid procedure satis- 
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fies detailed balance and therefore produces the correct 
statistical sampling 20] . 



III. APPLICATION TO PATH INTEGRAL 
SIMULATIONS 

In this section, we discuss the application of the 
method in Sect.[n]to path integral simulations. Here, we 
assume a 1-dimensional particle in a potential V. Gen- 
eralization to higher dimensions and many particles is 
straightforward. 



A. Definition of Length Scales 

Before describing the unigrid and multigrid algorithms, 
we define the concept of length scales in a path inte- 
gral simulation. Let the number of beads on the ring 
P be equal to 2 L , where L is a positive integer. (Since 
the path is cyclic, bead is identical to bead 2 L .) We 
divide the beads into different "levels" £ = 0, 1, ■■■L, 
such that I={lx2',3x2 £ ,5x 2 e , ■ • ■}. For examples, 
{1, 3, 5, ■ • ■} would belong to I = 0, {2, 6, 10, • • •} to I = 1, 
{4, 12,20, •••} to i = 2, etc. 

Using these definitions, we say the full path is of length 
scale L. Starting with the full path, we divide it into 
segments of different length scales at specific endpoints. 
Bisecting the full path yields two equal-length segments 
of length scale L — 1, the first one having endpoints 
and 2 L_1 and the second 2 L_1 and 2 L . Bisecting each 
of these L — 1 length scale segments again, we get four 
segments of length scale L — 2, and so on. 

In the absence of a potential V, every path segment of 
every length scale can be sampled independently from a 
gaussian distribution. Therefore, even though beads on 
different levels i are connected with each other via the 
kinetic energy springs, the path segments on all length 
scales can actually be generated in a completely uncor- 
related manner. But the presence of a nonzero potential 
V introduces correlations back into the path segments, 
and they can no longer be sampled independently. The 
potential produces a "confinement" effect on the path, 
which couples path segments of certain length scales. The 
length scale of these correlations depends on the spatial 
extent of the confining potential as well as the tempera- 
ture. For a fixed temperature, a broad and shallow po- 
tential produces less correlation than a steep narrow po- 
tential. More complicated potentials (those frequently 
present in condensed systems) may produce correlations 
on multiple length scales. For example, a bistable poten- 
tial produces two length scales, one for intra-well quan- 
tum fluctuations and the other for inter-well fluctuations. 
The goal of the stochastic blocking method is to remove 
some of these correlations produced by the potential, and 
the multigrid method furthermore refines it by attempt- 
ing to remove these correlations over all length scales. 



(a) (b) (c) (d) 
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FIG. 1: Illustration of the unigrid stochastic blocking proce- 
dure. (See text for details.) 



B. Stochastic Blocking: The Unigrid Method 

The stochastic blocking method aims at killing the cor- 
relations that come from the potential terms V(xj) in 
Eqn.JIJ. Following the ideas described in Sect.[n] we kill 
every V(xj) by either deleting it or freezing it. This is 
illustrated in Fig. ^ f° r a path with L — 5. Panel (a) 
shows the original path, with the V(xj) on all beads i 
depicted pictorially as open circles. The stochastic block- 
ing is done by first freezing the endpoints of the full path 
at j = = 2 L and then killing every potential term 
V(xj) for < j < 2 L . This operation is represented in 
panel (b) . The dotted lines highlight the beads on which 
the killing operation is performed, and the result of each 
killing operation is either a frozen bead (represented by 
a closed circle) or a deleted bead (represented by the 
absence of a circle). The path breaks up into frozen seg- 
ments consisting of beads 0(= 32), 3-6, 9, 16-23 and 
29-30 and the intervening deleted segments. The frozen 
segments can not be moved, but the deleted segments 
can be sampled independently from gaussian distribu- 
tions for free-particle paths of various lengths. Panel (c) 
shows the new path after the move (solid line) compared 
to the initial path (dashed line). All the potential terms 
V(xj) are finally restored resulting in the final path in 
panel (d). This completes one pass and the next pass 
begins anew with the killing procedure performed on the 
path in panel (d). (Since the origin of the cyclic path at 
j = = 2 L is always frozen in this method, we select 
a different origin at random on every pass to maintain 
ergodicity.) 

Before moving on to the multigrid implementation, 
we want to mention one minor difference between what 
is described here and the original implementation of 
the stochastic blocking procedure due to Swendson and 
Wang 01 as summarized in Sect. |nj The stochastic 
blocking procedure was originally applied to the Ising 
model, a system where the state of each particle belongs 
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FIG. 2: Illustration of the multigrid procedure. (See text for 
details.) 



to a finite discrete set and every potential term u a in 
the action of the Ising model is bounded from above by 
it*. As discussed in Sect. [HI to ensure that the deletion 
probabilities pd € [0,1], c a can be chosen to be less than 
cxp(— m* ) for every potential term. In our path integral 
application, however, the coordinate of each bead is a 
continuous variable and the potential is in general not 
bounded from above. This means that to strictly satisfy 
the requirement that pd < 1, c Q must be chosen to be 
which then results in no deletion at all. In practice, 
we can circumvent this minor problem by choosing a suf- 
ficiently large u* a and monitors the frequency at which 
the bound pd < 1 is violated during the simulation. By 
adjusting u* a to yield a bound violation frequency of less 
than 0.01%, we can maintain a relatively high deletion 
ratio while introducing negligible errors to the results. 



C. The Multigrid Method 

The multigrid method makes use of the stochastic 
blocking procedure in the unigrid method, but forces the 
path segments to break up on a predefined set of length 
scales. The procedure is illustrated pictorially for a path 
with L = 3 in Fig. 



Starting with the initial path in panel (a), we kill all 
V(xj) on levels I > 2. This results in a frozen bead at j — 
4. The open circles in panel (b) indicate V(xj) that are 
kept alive. When the new path segments are generated, 
they must be accepted or rejected based on a Metropolis 
criterion involving all the live beads. Panel (c) shows the 
new path after the move: the new segment between j = 
and 4 is accepted based on the three live beads at j = 1-3, 
but the other new segment between j = 4 and 8 is rejected 
based on the three live beads at j = 5-7; as a result, the 
new path (solid line) between j = 4 and 8 coincides with 
and the old one (dashed line). Panel (d) shows the new 
path on this level after all V{xj) are restored. Then we 
move on to the next finer level. On this level, we kill all 
beads on level I > 1, namely at j = 2, 4 and 6. This 
results in the frozen (solid circles), deleted (no circles) 
and live (open circles) beads indicated in panel (f). Two 
new segments are generated independently and accepted 
based on the live bonds at j = 1 and j = 3, 5 and 7. 
All beads are then restored and the procedure repeats on 
the finest level in the bottom row of Fig- El starting with 
killing all V{xj) for I > 0. 

Notice that the method presented here is very different 
from a method previously described by Janke and Sauer 
[HI ITi| , which they also call a "multigrid path integral 
method" . 

IV. RESULTS 

We have carried out numerical tests on our algorithm 
for several 1-dimensional quantum systems and compared 
them against other methods, including: 

1. Metropolis: a conventional Metropolis algorithm 
based on single-bead moves in the ir-coordinates; 

2. Bisection: the bisection algorithm as described by 
Ceperley in Q; 

3. Unigrid: the unigrid algorithm as described above; 
and 

4. The so-called "multi grid path integral method" of 
Janke and Sauer (l5lll6| . 

The tests were carried out on different model systems 
with symmetric single- and double- well potentials. In the 
simulations, we have adopted a system of dimensionless 
units in which m = Ti = 1. For all tests, (3 — 10, yielding 
a thermal wavelength of approximately 1.6 which sets the 
length scale of the quantum dispersion due to the kinetic 
part of the action. For all of the systems studied, the 
ground state dominates at this temperature. Because the 
potentials are all symmetric, (x) should be 0; therefore, 
how quickly the measured {x) goes to the exact value of 
will provide a good estimate of the efficiency of each 
algorithm. Alternatively, if we calculate the quantity x = 
Y^j=i x j after each Monte Carlo step, we can determine 
the efficiency of the algorithm by examining how rapidly 
x fluctuates around the exact value (x) = 0. 
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FIG. 3: Initial and final paths after 10 4 MCS in a conven- 
tional Metropolis Monte Carlo simulation of Model A using a 
discretization of L = 12. The configuration has hardly moved, 
indicating the severity of the slowing-down problem. 
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A. Harmonic Potential (Model A): V{x) = \x 2 

Model A is a simple harmonic potential. Figure |3| 
shows the initial and final quantum paths after 10 4 single- 
particle Metropolis Monte Carlo steps (one MCS is de- 
fined as having every bead on the path subjected to one 
trial move on the average). There is very little move- 
ment in the configuration of the whole path. Clearly, 
the conventional Metropolis algorithm is highly ineffec- 
tive. P — 4096 beads, or L = 12, were used to represent 
the path here. With P of this magnitude, the bead-to- 
bead dispersion is merely 0.05, roughly 3% of the thermal 
wavelength, making the harmonic bonds between succes- 
sive beads along the ring extremely stiff. 

The discretization L = 12 is much larger than what 
is needed for the path integral results to converge to the 
continuum limit for this model. The minimum required 
L is 5. Figure 0] shows i as a function of MCS for the 
four methods for L — 5. Visually, we can see that the 
multigrid method is the most efficient, the unigrid and 
the bisection methods are comparable and slower than 
the multigrid method, and the Metropolis method is the 
least efficient. To get a more precise measure of the effi- 
ciencies, the data in Fig. 0] were autocorrelated and the 
correlation functions are shown in Fig. [S] for the four 
methods. The general conclusions we obtained from a 
visual inspection of Fig. 0] are confirmed by Fig- El — the 
multigrid method has the fastest decay time constant, 
approximately 1 MCS, and therefore is the most efficient. 

The efficiencies of the four different methods for other 
values of L are shown in Fig. The left panel of Fig. 
shows the decay time constants t c in units of MCS from 
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FIG. 5: Autocorrelation functions of x for the data in Fig. [I] 
The multigrid method has the fastest decay time. 



the autocorrelation function of x for the four methods at 
different path discretizations L. The vertical line at L — 
5 indicates the minimum value of L needed for the path 
integral results to converge for this model. Everything 
on the left of this line does not converge to the correct 
continuum limit. Therefore, the only relevant results are 
those to the right of this line. 

From a computational standpoint, the ultimate mea- 
sure of real efficiency is the correlation time measured in 
actual CPU cycles (i.e. r c in units of MCS multiplied 
by CPU time per MCS) since difference algorithms will 
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FIG. 6: Correlation times r c from four different MC meth- 
ods at different discretizations L for Model A. The vertical 
line at L — 5 indicates the minimum L needed for the path 
integral results to converge to the continuum limit. The left 
panel shows correlation times in units of MCS. The right panel 
shows the same but in units of actual CPU millisecond. (The 
segment length used for the bisection method was 2 L ~ 2 for 
each L.) 



incur different CPU time per MCS. The right panel in 
Fig. H3 shows correlation times in actual CPU millisec- 
onds. The real efficiency of the multigrid method at the 
minimum required discretization for convergence (L = 5) 
is about a factor of 4 better than its nearest competitor, 
the bisection method, and a factor of 40 better than the 
standard Metropolis method. 

This model has essentially one length scale dictated by 
the confinement effect of the potential. For the particular 
parameters chosen for this model, the confinement length 
scale of the potential turns out to be quite similar to the 
natural thermal wavelength of the free-particle path. For 
this situation, a relatively small L is sufficient for the path 
integral to converge to its continuum limit. Moreover, 
the four different methods, with their correlation times 
spanning only one and a half orders of magnitude, do not 
exhibit vastly different efficiencies. 

These conclusions drawn from Model A are certainly 
not universal. The relative efficiencies of the various 
methods will depend on a number of factors, such as the 
length scale of the potential and the temperature. The 
additional examples presented below demonstrate this. 
But in all the models considered, the multigrid method 
was always the most efficient by an order of magnitude 
or more compared to any other method. 



B. Compressed Harmonic Potential (Model B): 

V(x) = \{ixf 

In Model B, the length scale of the potential is sub- 
stantially smaller than the thermal wavelength. This sit- 



FIG. 7: Correlation times r c from four different MC methods 
at different discretizations L for Model B. The minimum dis- 
cretization required for convergence is L = 7, indicated by the 
vertical line. Left and right panels show r c in units of MCS 
and CPU millisecond, respectively. (The segment length used 
for the bisection method was 2 L ~ 3 for each L.) 



uation is quite typical in real condensed-phase quantum 
systems. The decay time constants of the autocorrelation 
functions of x are summarized in Fig.0for different L and 
the four MC methods. For this model, L = 7 is the mini- 
mum required for convergent path integral results. At or 
above this value of L, the multigrid method outperforms 
all the other methods in both MCS and CPU efficiencies. 
For this model, the second most efficient method is the 
unigrid method, which performs at about a factor of 4 
poorer than the multigrid method at L = 7 in terms of 
real CPU efficiency. On the other hand, the Metropolis 
method is a factor of 80 less efficient. 



C. Double- Well Potential (Model C): 

V(x) = ~3x 2 + x 4 

Model C is a double-well potential with a moderate 
barrier. The inter-well separation is somewhat longer 
than the thermal wavelength but not by much. This is 
the first model that has at least two length scales due 
to the bistable nature of the potential. Because of the 
moderate barrier, the length scales of the intra-well and 
inter-well quantum fluctuations are in reality not very 
different from each other. 

The decay time constants of the autocorrelation func- 
tions of x are summarized in Fig. [S] for different L and 
the four MC methods. For this model, L = 7 is the mini- 
mum required for convergent path integral results. At or 
above this value of L, the multigrid method outperforms 
all the other methods in both MCS and CPU efficiencies. 
For this model, the next most efficient method is the bi- 
section method, which performs at about a factor of 3 
poorer than the multigrid method at L = 7 in terms of 
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FIG. 8: Correlation times r c from four different MC methods 
at different discretizations L for Model C. The minimum dis- 
cretization required for convergence is L — 7, indicated by the 
vertical line. Left and right panels show r c in units of MCS 
and CPU millisecond, respectively. (The segment length used 
for the bisection method was 2 L ~ 2 for each L.) 



real CPU efficiency. 



D. Compressed Double-Well Potential (Model D): 

V(x) = -3(2x) 2 + (2a;) 4 

Model D has the same double-well potential as 
Model C, but the potential is compressed in the IE- 
direction. In addition to having two distinct length 
scales, this model is further complicated by having a 
rather severe confinement effect because the length scale 
of the potential is much smaller than the natural thermal 
wavelength of the free-particle path. The result is that 
the minimum discretization required for convergence be- 
comes larger [L = 8) and the multigrid method becomes 
increasingly advantageous compared to the other meth- 
ods. On the other hand, the bisection method suffers 
here, because it generates path segment of only one pre- 
defined length scale. When that length scale is adjusted 
for maximum optimal efficiency, it would match one but 
misses all the other relevant length scales present in the 
problem plj . 

The decay time constants of the autocorrelation func- 
tions of x are summarized in Fig. [5] for different L and 
three MC methods. The results from the unigrid method 
are not shown because for every L the decay time in the 
unigrid results is greater than 10 MCS and we were un- 
able to equilibrate the unigrid simulations. The multigrid 
method outperforms the other two methods in both MCS 
and CPU efficiencies. For this model, the next most effi- 
cient method is the bisection method, which performs at 
about a factor of 10 poorer than the multigrid method 
at L = 8 in terms of real CPU efficiency. 
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FIG. 9: Correlation times r c from three different MC meth- 
ods at different discretizations L for Model D. The minimum 
discretization required for convergence is L = 8, indicated by 
the vertical line. Left and right panels show r c in units of MCS 
and CPU millisecond, respectively. (The segment length used 
for the bisection method was 2 L ~ 4 for each L.) 



E. Model of Janke and Sauer (Model E): 

V(x) = -0.5x 2 + O.Oix 4 

To make contact with the results of Janke and Sauer 
|l5l ITtH , we carried out simulations on the double- well 
potential studied in their papers, using the same param- 
eters they have used. The decay time constants of the 
autocorrelation functions of x are summarized in Fig. 1101 
for different L. L = 7 is the minimum required for con- 
vergent path integral results. The Metropolis, unigrid 
and bisection results are not shown because their decay 
times are all greater than 10 4 MCS for all values of L. 
At L = 7, the multigrid method is about an order of 
magnitude more efficient than the Janke/Sauer method. 
But notice that the CPU correlation time for this model 
is much larger (r c ss 100 CPU ms) compared to all the 
previous models. Therefore, this model potential seems 
to present a slightly more challenging problem even for 
the multigrid method. 



F. An Electron in the Field of Two Positive Ions 
(Model F) 

This last model consists of a 1-dimensional electron in 
the field of two ions of +2 charge separated by a dis- 
tance of 15 A. To prevent the electron path from collaps- 
ing onto either ion, we set up a 0.25 A hard core radius 
around each. The temperature is 300 K and the dielec- 
tric constant is 78. Similar to some of the other models 
already considered, this model consists of a bistable po- 
tential with a high barrier between the two wells. This 
model is however qualitative quite different from the oth- 
ers, because the coulombic potential centered on each ion 
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FIG. 10: Correlation times r c from three different MC meth- 
ods at different discretizations L for Model E. The minimum 
discretization required for convergence is L — 7, indicated by 
the vertical line. Left and right panels show r c in units of 
MCS and CPU millisecond, resnectivelv. 
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FIG. 11: Measurement of x after each MCS during the course 
of the MC simulation for Model F using the Metropolis, the 
bisection, the unigrid and the multigrid methods using a dis- 
cretization of L — 9. 



is much narrower in comparison with the fairly open well 
bottoms in the other models. This results in a strong 
confinement effect on the electron path. The natural 
quantum fluctuations of the path are of a longer scale 
than the width of the potential wells. Figure ITT1 shows x 
as a function of MCS for the four methods at L = 9, the 
minimum required discretization needed for convergence. 



Both the bisection and the Metropolis methods failed to 
equilibrium in 10000 MCS. The decay time constants of 
the autocorrelation function of x are shown in Fig.ll2lfor 

the : 
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FIG. 12: Correlation times t c from the multigrid and the 
unigrid methods for Model F. The minimum discretization 
required for convergence is L — 9, indicated by the vertical 
line. Left and right panels show r c in units of MCS and CPU 
millisecond, respectively. 



V. CONCLUSION 

We have described a new Monte Carlo method for sim- 
ulating imaginary-time path integrals. The new method 
uses a combination of a stochastic blocking algorithm and 
multigrid ideas to completely eliminate the slowing-down 
problem in the sampling of discretized quantum paths 
near the continuum limit. The method has been tested 
on several 1-dimensional quantum systems and found to 
exhibit highly ergodic dynamics. The new method offers 
distinct advantages over other methods in cases where the 
length scale of the potential is smaller than the quantum 
dispersion of the path, a situation that is quite typical 
in real condensed-phase quantum systems. On the other 
hand, for systems with bistable potentials with length 
scales much larger than the quantum dispersion of the 
path, the new method, though better than all the other 
methods, has only limited utility. 
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